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Abstract 

Until  recently,  first-principles  calculations  of 
potential  energy’  surfaces  ( PES )  were  restricted  to 
intermolecular  interactions  involving  molecules 
containing  just  several  atoms.  On  one  hand,  this  was  due 
to  the  high  cost  of  wave-function-based  electronic 
structure  methods  and,  on  the  other  hand,  to  the  failure  of 
the  density  functional  theory  (DFT)  approaches  to 
reproduce  the  dispersion  part  of  intermolecular 
interactions. 

One  solution  to  this  problem  is  symmetry-adapted 
perturbation  theory’  based  on  DFT  description  of 
monomers  [SAPT(DFT)].  In  applications  to  energetic 
materials,  SAPT(DFT)  predicted  the  correct  crystal 
structure  of  RDX  (l,3,5^trinitroperhydro-l,3,5-triazine). 
Recently,  the  complete  PES  of  FOX-7  (1,1 -diamino-2, 2- 
dinitroethene)  dimer  was  obtained  using  SAPT(DFT). 
Preliminary’  molecular  dynamics  simulations  of  the  FOX- 
7  crystal  show  an  improved  agreement  with  experiment 
compared  to  literature  results.  A  recently  developed 
nearly-linear  scaling  implementation  of  the  SAPT(DFT) 
dispersion  energy’  has  been  applied  to  interactions  of 
energetic  molecules.  When  the  development  of  linear- 
scaling  SAPT(DFT)  is  finished,  accurate  studies  of 
energetic  molecules  significantly  larger  than  RDX  and  of 
other  important  systems  (including  biomolecules), 
containing  in  excess  of  one  hundred  atoms,  will  be 
possible. 

Another  approach  which  can  be  applied  to  such 
systems  is  the  dispersion-less  density  functional  (dlDF) 
method  developed  in  our  group  which  reproduces 
interaction  energies  with  the  dispersion  component 
removed.  The  dispersion  energy’  is  then  computed  from 
an  asymptotic  function  fitted  to  SAPT(DFT)  dispersion 


energies  of  a  training  set,  resulting  in  a  method  denoted 
as  dlDF+D.  Cross  sections  of  the  PES  of  the  HMX 
(octahydro-1,3,5, 7-tetranitro-l,3,5, 7-tetrazocane)  dimer 
calculated  using  dlDF+D  are  presented. 

1,  Introduction 

Energetic  materials  are  important  to  the  national 
defense,  and  the  study  of  their  bulk  properties  represents 
one  of  the  most  active  fields  of  the  material  sciences. 
Some  of  the  obvious  disadvantages  of  experimental 
research  in  this  area  are  the  relatively  high  risk  involved 
and  associated  costs,  especially  when  synthesizing  and 
testing  new  materials  of  yet  unknown  properties. 
Therefore,  a  reliable  theoretical  prediction  of  the  crystal 
structure  and  other  physicochemical  properties  of  such 
compounds  is  highly  desirable.  The  most  important 
prerequisite  for  achieving  this  goal  is  the  availability  of  an 
intermolecular  potential  (force  field)  representing  the 
interaction  energy  of  a  pair  of  molecules  as  a  function  of 
their  separation  and  orientation.  One  possible  way  to 
develop  such  potentials  is  to  fit  their  functional  form 
(including  a  number  of  adjustable  parameters)  so  as  to 
reproduce  as  closely  as  possible  certain  experimental  data. 
The  most  serious  drawback  of  empirical  force  fields  is 
their  uncertain  predictive  power.  Since  such  fields  only 
contain  information  about  the  pertinent  system  under  the 
conditions  corresponding  to  the  data  used  in  the  fitting, 
these  potentials  often  fail  when  applied  to  significantly 
different  conditions  or  to  properties  not  related  to  those 
data.  On  the  other  hand,  ab  initio  potentials  are  obtained 
from  first-principles  quantum  mechanical  calculations  and 
can  be,  formally,  as  accurate  as  desired,  provided  that  a 
sufficiently  large  basis  set  and  a  high  level  of  theory, 


especially  an  appropriate  description  of  the  electron 
correlation,  are  used.  Wave-function-based  methods, 
such  as  the  coupled  cluster  theory  or  symmetry-adapted 
perturbation  theory  (SAPT),  are  quite  demanding 
computationally,  and  their  cost  scales  too  steeply  with  the 
system  size  to  be  applicable  to  accurate  calculations  of 
molecules  larger  than  about  a  dozen  of  atoms.  Density 
functional  theory  (DFT),  while  significantly  less 
expensive,  poorly  describes  long-range  dynamic 
correlation  effects,  and  thus  the  dispersion  interaction, 
which  is  often  the  predominant  effect  in  compounds 
forming  molecular  crystals.  Because  of  this  dilemma,  a 
successful  first-principles  simulation  of  energetic 
materials  was  not  possible  as  recently  as  a  decade  ago. 
The  situation  changed  with  an  advent  of  new  methods 
combining  the  advantages  of  various  distinct  theories. 
The  symmetry-adapted  perturbation  theory  with  the 
density-functional  theory  description  of  monomers11-71, 
abbreviated  as  SAPT(DFT)  or  DFT-SAPT,  is  similar  to 
the  regular  SAPT  approach,  but  represents  the  monomers 
using  the  Kohn-Sham  DFT  approach1*1  instead  of  the 
much  more  expensive  coupled  cluster  theory.  The  crucial 
dispersion  interaction  is  calculated  from  the  so-called 
frequency-dependent  density  susceptibilities  (FDDS)  (i.e., 
polarization  propagators)  of  the  monomers14'51,  which  in 
turn  can  be  obtained  from  a  coupled  Kohn-Sham  (time- 
dependent  DFT)  calculation.  A  big  advantage  of  the 
SAPT(DFT)  approach  is  the  fact  that  the  full  dimer 
calculation  is  entirely  avoided  and  the  once-computed 
FDDS’s  of  the  monomers  can  be  used  for  any  dimer 
configuration,  provided  that  the  monomers  are  rigid.  The 
density-fitting  technique[9]  has  been  implemented  in  this 
method1  IH  15  and  has  resulted  in  an  additional 
improvement  of  the  efficiency.  The  SAPT(DFT) 
approach  was  applied  to  such  computationally  demanding 
problems  as  the  complete  PES  of  the  benzene  dimer[13] 
and  selected  configurations  of  the  pyrene  dimer[  ,  the 
latter  containing  as  many  as  52  atoms.  In  the  field  of 
energetic  materials,  SAPT(DFT)  has  been  very 
successful115-171  in  predicting  the  correct  crystal  structure 
of  RDX  (l,3,5-trinitroperhydro-l,3,5-triazine),  with  the 
dimer  of  RDX  containing  42  atoms.  We  also  calculated 
the  complete  PES  of  the  FOX-7  (1,1 -diamino-2, 2- 
dinitroethene)  dimer  and  here  we  report  some  initial 
results  of  the  simulations  for  the  FOX-7  crystal. 

The  application  range  of  SAPT(DFT)  can  be 
extended  to  still  larger  dimers,  of  the  size  in  excess  of  a 
hundred  of  atoms,  by  the  near-linear  scaling  version,  L- 
SAPT(DFT).  In  this  approach,  one  divides  the  set  of 
atoms  in  each  monomer  into  a  “far-region”  and  a  “near¬ 
region”  subsets.  The  far-region  atoms  are  those  whose 
distance  from  every  atom  in  the  other  monomer  is  larger 
than  a  certain  cutoff  value,  i?cut.  Obviously,  for  very  large 
molecules  most  atoms  belong  to  the  far  region.  Since  the 
dimer  interaction  can  be  represented  as  a  sum  of 


interactions  between  disjoint  regions  of  both  monomers, 
only  a  small  fraction  of  the  total  calculation  (between 
near-region  atoms  of  both  monomers)  has  to  be  performed 
using  the  regular  SAPT(DFT)  expressions,  while  the  rest 
can  be  sufficiently  accurately  evaluated  from  simple 
asymptotic  formulas.  The  results  of  L-SAPT(DFT) 
calculations  have  already  been  reported  for  the 
electrostatic  energy1181  and  the  strategy  described  above 
has  recently  been  implemented  for  the  dispersion  energy. 
Our  initial  results  for  the  latter  component  are  presented 
below. 

Recently,  our  group  introduced1191  another  new 
hybrid  method  for  calculations  of  interaction  energies 
with  a  large  dispersion  component.  The  idea  of  this 
method  is  based  on  the  observation  that  DFT  should  be 
able  to  quite  accurately  predict  all  the  components  of  the 
interaction  energies  but  the  dispersion  energy. 
Consequently,  a  “dispersionless”  density  functional 
(dlDF)  has  been  optimized1191  to  reproduce  the 
dispersionless  interaction  energies,  i.e.,  the  interaction 
energies  with  the  dispersion  components  subtracted,  for  a 
training  set  containing  nine  common  dimers  at  various 
intermonomer  separations.  The  interaction  energies  were 
obtained  using  the  coupled  cluster  method  with  single, 
double,  and  noniterated  triple  excitations  [CCSD(T)], 
while  the  dispersion  energies  were  represented  by  the  sum 
of  the  second-order  dispersion  and  exchange  dispersion 
energies,  E17,1  +  CcL&p  >  calculated  with  SAPT (DFT). 

The  interaction  energies  produced  by  dlDF  calculations 
must,  of  course,  be  supplemented  with  dispersion  energies 
obtained  by  some  other  method,  forming  a  hybrid 
approach  known  as  dlDFT+D.  Up  to  some  dimer  size, 
one  can  use  the  dispersion  energies  computed  using 
SAPT(DFT).  For  larger  systems,  the  authors  of 
Reference  19  developed  a  simple  damped  asymptotic 
dispersion  function,  in  the  form  of  a  sum  of  pair  (atom- 
atom)  contributions.  Its  parameters  were  optimized  to 
reproduce  the  sum  +  E^h_disp  on  a  training  set.  The 

dispersion  function  was  significantly  improved  in 
Reference  20  by  introducing  different  parameter  values 
for  hydrogen  atoms,  depending  on  their  connectivity,  and 
by  enlarging  the  training  set.  This  form  of  the  dlDFT+D 
approach  was  used  in  the  present  work. 

2.  FOX-7  Calculations 

The  SAPT(DFT)  method  was  recently  used  to 
calculate  about  1,300  points  of  the  FOX-7  dimer  PES.  A 
pair  potential  based  on  Coulomb  terms  describing  the 
interaction  between  partial  charges  on  atoms  plus 
Buckingham-type  exp-6  terms  was  obtained.  The  partial 
charges  were  determined  by  fitting  to  ab  initio  molecular 
multipole  moments  of  the  monomer  (through  1=6).  The 


remaining  parameters  were  fitted  to  the  SAPT(DFT) 
interaction  energies  directly.  Here  we  present  (see  Table 
1),  for  the  first  time,  results  of  a  NsT  (constant  number  of 
molecules,  stress,  and  temperature)  simulation  of  the 
crystal  lattice  parameters,  based  on  the  new  potential, 
compared  to  the  semi-empirical  potential  of  Sorescu,  et 
al.[:l].  Except  for  the  lattice  dimension  a,  the  SAPT(DFT) 
potential  clearly  improves  the  agreement  with 
experimental  data  and  the  error  of  the  unit  cell  volume  is 
reduced  by  a  factor  of  1.7.  The  work  on  FOX-7  is  in 
progress.  In  particular,  calculating  new  PES  points  in 
some  critical  regions  is  still  an  option,  which  would  likely 
further  improve  the  agreement  with  experiment. 


Table  1.  Results  of  the  NsT  simulation  of  the  FOX-7 
crystal  lattice  at  300K.  The  unit  cell  parameters  are  in 
angstroms  (a,  b,  c)  and  degrees  (a,  p,  j),  the  volume  in 
A\ 


method  a  b  c  a  f)  7  volume  (%error) 

Experiment  [22|  6.93  6.62  11.31  90.00  90.06  90.00  5l9 

Ref.  [21)  6.93  6.86  11.60  89.99  89.79  89.98  551  (6.2) 

SAPT(DFT)  6,79  6.54  1 1.25  89,99  90.02  89.98  500  (3.7) 


3.  HMX  Dimer  Calculations 

The  HMX  (octahydro-1, 3,5, 7-tetranitro-l, 3,5,7- 
tetrazocane)  molecule  is  an  eight-membered  analogue  of 
the  six-membered  RDX  ring  and  its  dimer  has  56  atoms. 
HMX  exists  in  four  different  polymorphic  forms,  labeled 
a  through  5.  The  form  stable  at  ambient  conditions,  //- 
HMX,  crystallizes  in  the  monoclinic  space  group  P2i/c.  It 
is  the  only  crystal  structure  of  HMX  containing  molecules 
in  the  Crsymmetry  conformation  (the  three  other  forms 
contain  CVsymmetric  molecules).  The  bulk  properties  of 
HMX  crystals  were  the  subject  of  numerous  atomistic 
simulation  papers[2V3S1.  However,  to  our  knowledge, 
there  is  no  accurate  study  of  the  dimer  in  the  literature. 
We  performed  such  a  study  applying  the  dlDF+D  method. 
First,  we  used  the  simple  HMX  pair  potential  developed 
in  Reference  24  to  perform  a  quick  scan  of  the  PES.  The 
P- HMX  monomer  geometry  was  retrieved  from  the 
Cambridge  Structural  Database  (code  OCHTET12).  It 
was  kept  frozen  throughout  all  the  calculations,  i.e.,  only 
the  rigid-monomer  dimer  configurations  were  studied. 
We  calculated  the  dimer  interaction  energies  on  the  six¬ 
dimensional  grid  constructed  by  taking  all  the  center-of- 
mass  (COM)  distances  R  from  5  to  7.25  A  with  a  step  of 
0.25  A  and  all  the  possible  relative  Euler  angles  with  a 
step  of  10°.  The  geometry  corresponding  to  the  lowest 
energy  was  subsequently  refined  in  several  steps  by 
forming  denser  and  denser  grids  around  the  current 
minimum.  The  obtained  minimum  energy  structure 
(interaction  energy  equal  to  -5.86  kcal/mol)  is  presented 
in  Figure  1.  It  has  the  COM  distance  of  5.42  A  and  both 


molecules  are  almost  exactly  parallel  (with  a  small 
deviation  of  4°  caused  probably  by  the  fact  that  the 
monomer  of  Reference  24  does  not  have  the  exact  C, 
symmetry). 


Figure  1.  Minimum  energy  geometry  of  the  HMX  dimer 


In  the  next  step,  we  used  the  dlDFT+D  method  and 
the  correlation-consistent  augmented  double-zeta 
Dunning  basis  set  (aug-cc-pVDZ,  1 064  basis  functions)1291 
to  obtain  a  cross  section  of  the  HMX  dimer  PES, 
containing  the  structure  of  Figure  1  and  thirteen  other 
structures  with  the  same  relative  orientation  and  with  R 
ranging  from  4.75  to  8  A.  The  summary  of  the  results  is 
presented  in  Table  2.  Qualitatively,  the  shape  of  both 
interaction  potentials,  empirical  and  dlDF+D,  along  this 
cross  section  is  similar.  However,  the  dlDFT+D 
calculation  predicts  a  significantly  deeper  minimum 
(-7.00  vs.  -5.86  kcal/mol)  at  a  slightly  shorter  COM 
distance  (5.30  vs.  5.42  A).  The  interaction  is  strongly 
dominated  by  the  dispersion  energy,  which  is  (at  the 
minimum  distance)  2.5  times  larger  in  absolute  value  than 
the  sum  of  the  remaining  components. 


Table  2.  HMX  dimer  interaction  energies  Emt  (in 
kcal/mol)  for  the  orientation  of  the  Sewell’s 
potential1241  global  minimum  as  functions  of  the 
center-of-mass  distance  between  the  monomers,  R(in 
angstroms).  Erep  and  Eeiec  stand  for  the  repulsion  and 
electrostatic  energy  terms  in  the  potential  of 
Reference  24  and  EdiSp  for  the  dispersion  energy 
calculated  with  either  method. 


R 

dlDFT+D  method 
dIDFT  C,h,„  £),» 

potential  of  Ref.  |24] 

Erep  +  Ev\vc  Edlsp  E(nt 

4.75 

24.1528 

-24.7671  -0.6143 

29.0331 

-24.5247  4.5084 

5.00 

1  1.0230 

-17.4288  -5.5058 

14.3095 

-17.2939  -2.9844 

5.20 

6.4220 

-13.2496  -6.8276 

7.9889 

-13.2431  -5.2543 

5.25 

5.4408 

-12.3852  -6.9444 

6.8905 

-12.4098  -5.5193 

5.30 

1.5832 

-11.5824  -6.9993 

5.9382 

-11.6369  -5.6987 

5.35 

3.8384 

-10.8367  -6.9983 

5.1133 

-10.9194  -  5.8061 

5.42 

2.9629 

-9.8804  -6.9176 

4.1417 

-9.9998  -5.8581 

5.50 

2.1521 

-8.9008  -6.7487 

3.2492 

-9.0575  -5.8083 

5.70 

0.7586 

-6.8945  -6.1359 

1.7567 

-7.1227  -5.3660 

6.00 

-0.2075 

-4.7742  -4.9817 

0.6843 

-5.0573  -4.3729 

6.50 

-0.6952 

-2.6995  -3.3947 

0.1356 

-2.9839  -2.8482 

7.00 

-0.7372 

-1.6057  -2.3429 

0.0254 

—  1.8442  -1.8187 

7.50 

-0.6655 

-0.9993  -1.6648 

0.0044 

-1.1850  -1.1806 

8.00 

-0.5764 

-0.6467  -1.2231 

0.0005 

-0.7869  -0.7864 

Table  3  shows  results  of  the  application  of  the  new, 
near-linear  scaling  dispersion  code  of  the  SAPT(DFT) 
method  to  one  of  the  configurations  from  Table  2  (R= 5.42 
A).  The  basis  set  of  Sadlej[40]  was  used,  containing  552 
basis  functions  per  monomer.  The  calculation  with  the 
largest  cutoff  distance,  Rcm=6.()  A,  assigns  all  the  atoms  to 
the  near  region,  which  results  in  the  application  of  the 
regular  SAPT(DFT)  algorithm  and  yields  the  exact 
(within  this  particular  basis  set)  value  of  the  dispersion 
energy,  Aii-,P=~8.88  kcal/mol.  This  benchmark  result  can 
be  compared  with  the  value  from  the  damped  asymptotic 
dispersion  function  of  the  dlDFT+D  method,  equal  to 
-9.88  kcal/mol  (11%  of  error).  In  fact,  the  true 
discrepancy  is  somewhat  larger  since  the  former  value 

does  not  include  -E^cVdisp ,  which  is  positive.  Reducing 

the  cutoff  distance  introduces  a  gradually  increasing  error 
of  the  L-SAPT(DFT)  dispersion  energy,  but  up  to  Rcut=4 
A,  where  half  of  all  the  atoms  are  in  the  far  region,  L- 
SAPT(DFT)  predicts  the  dispersion  energy  more 
accurately  than  the  current  version[2n]  of  the  dlDFT+D 
dispersion  fit.  It  should  be  noted  that  this  particular 
system  represents  a  rather  unfavorable  case  for  a  L- 
SAPT(DFT)  calculation.  As  seen  in  Figure  1,  both 
molecules  are  stacked  in  a  near-sandwich  configuration, 
where  the  near  region  necessarily  extends  to  a  significant 
portion  of  both  molecules.  With  the  current  code,  a  single 
calculation  of  the  HMX  dimer  interaction  energy  requires 
between  2  and  3  CPU-days  with  either  dlDFT+D  or  L- 
SAPT(DFT)  on  the  MJM  linux  cluster  (3.0  GFlz  Intel 
Woodcrest  processors)  at  US  Army  Research  Laboratory 
(ARL). 

Table  3.  L-SAPT(DFT)  HMX  dimer  dispersion  energy 
Edisp  and  its  four  parts  (in  kcal/mol)  for  the  geometry  of 
the  Sewell’s  potential1241  global  minimum  as  a  function 

of  the  cutoff  distance,  RCut  (in  angstroms).  E^p 
denotes  the  contribution  from  the  near  region  (X=N)  or 
far  region  (X=F)  of  the  monomers.  E^J  was 

calculated  from  the  exact  SAPT(DFT)  formulas  and  the 
other  three  contributions  from  asymptotic  formulas. 
(N/N)  stands  for  the  number  of  near-region  atoms  in 
both  monomers. 

R  m  E™,  E^Z,  error  (%)  (N/N) 

3.0  -1.4988  -3.9151  -3.1285  -3.9320  -12.4744  -  40.5  (5/5) 

3.5  -3.6507  -3.4703  -2.6240  -1.3264  -11.0715  -24.7  (10/9) 

4.0  -6.6399  -2.1528  -0.5697  -0.2767  -9.6392  -8.6  (15/13) 

4.5  -8.5077  -0.1568  -0.2277  -0.0760  -8.9682  -1.05  (19/19) 

5.0  -8.8748  -0.0760  0.0287  -0.0393  -8.9614  -0.97  (23/22) 

5.5  -8.8774  -0.0037  0.0000  0.0000  -8.8811  -0.063  (28/26) 

6.0  -8.8754  0.0000  0.0000  0.0000  -8.8754  0.000  (28/28) 


4.  Conclusions 

The  recent  progress  in  developing  computationally 
efficient  first-principles  methods  of  studying  molecular 
interactions,  especially  in  dealing  with  the  dispersion 
energy  treatment,  has  made  possible  a  systematic  study  of 
dimers  where  first-principles  calculations  were 
prohibitively  time-consuming  just  a  few  years  ago.  In 
particular,  our  pilot  calculations  for  the  HMX  dimer  (56 
atoms)  show  that  obtaining  the  complete  PES  for  this 
important  energetic  material  is  now  feasible.  Assuming 
that  a  similar  number  of  PES  points  as  that  for  FOX-7  are 
required,  the  complete  project  will  consume  around  100K 
CPU-hours. 
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